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ABSTRACT 

We present a new, completely Lagrangian magnetohydrodynamics code that is based on the 
SPH method. The equations of self-gravitating hydrodynamics are derived self-consistently 
from a Lagrangian and account for variable smoothing length ("grad-h"-) terms in both the 
hydrodynamic and the gravitational acceleration equations. The evolution of the magnetic 
field is formulated in terms of so-called Euler potentials which are advected with the fluid and 
thus guarantee the MHD flux-freezing condition. This formulation is equivalent to a vector 
potential approach and therefore fulfills the V - B = O-constraint by construction. Extensive 
tests in one, two and three dimensions are presented. The tests demonstrate the excellent 
conservation properties of the code and show the clear superiority of the Euler potentials over 
earlier magnetic SPH formulations. 
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1 INTRODUCTION 

In many areas of astrophysics a transition from pure hy- 
drodynamic to magnetohydrodynamic simulations is underway. 
Magnetohydrodynamic calculations have a long-standing tradi- 
tion i n the context of c ore-co l lapse supernovae, see for ex - 
ample iLeBlanc & Wilson! Jl970h : iBisnovatvi-Kogan et al.l i fl976h : 
iMeier etalJ jl976t) : ISvmbalistvl dl984l) . In recent years MHD 
simu lations of core-collapse supernovae have seen a renaissance 
(e.g. lA kiyama et al. 2003; Mizuno et alj|2004l: iL iebendorfe r et al. 
2004 ; lYamada & Sawaill2004 iKotake et alj|2004 lArdelian etafl 
20051: IProga 2005; Ob ergaulinger et al . 2006; M asada et alj|2006l ; 
Shibat a et alj|2006l : iBurrows et alj|2007l) . mainly due to the con- 



clusion that a small fraction of core-collapse supernovae, those 
related to long gamma-ray bursts, require relativistic and well- 
collimated jets and due to the difficulty to make supernovae ex- 
plode via the delayed, neutrino-driven mechan ism. In accretion 
physi cs the magnetorotational instability (e.g. iBalbus & Hawleyl 
199a) is now the widely accepted mechanism behind the angular 
momentum transport that determines the accretion rate. Many of 
the recent accretion simulations are performed in the framework 
of (sometimes genera l relativistic) magnetohydrodynamics (e 
Stone & Prin"giefl200ll; Ibe Villiers et ail 12003; Isan'o et al.l [200- 



McKinnevI 120051 : lHawlev & Krolild [2006; McKinne v~& Naravanl 
20071) . Several of the more recent star and planet formation 
calcu lations have also included e ffects of the magnetic field 
(e.g. lHosking & WhitwortH |2004|; iNelsonl 120051: Izieglen 120051; 
Machi da et alj200d:lFromang & Nelsonll2006l ; lBaneriee & Pudritj 
20061 : IPrice&BatdlToOVbh . Most recently. MHD simulations 



were also used in t he context of com p act binary mergers (e.g . 
Shibata et aD 120061; | Duez et iD l2006l: IPrice & Rosswogl 120061 : 



iDuez et alj2007l) . 

Nearly all of the above calculations have been carried out on Eu- 
lerian gri ds. Some SPH-formulations that include magnetic fields 
exist (e.g. Gingold & Monaghanll977l:|Phillips & Monaghai j 19851 : 
lB0rve et alj200ll ; lDolag et al.l2002| ; |Price & Monaghanil2005l) . but 
obtaining a stable formulation has proved notoriously difficult, not 
least because of the difficulty in f ulfilling the V ■ B = 0-constraint 
on Lagrangian particles (see e.g. IPrice & Bate! |2007a) for a brief 
review). Nevertheless, for applications that involve large deforma- 
tions Lagrangian schemes have definite advantages. 
Here we present a detailed description of our Lagrangian magneto- 
hydrodynamics code, MAGMA ("a magnetohydrodynamics code 
for merger applications"). Our developments are mainly driven by 
the application to mergers of magnetized neutron stars, but the 
described methods are applicable to smoothed particle (magneto- 
) hydrodynamics in general. Ingredients of the code that have 
been described elsewhere are briefly summarized, the interested 
reader is referred to the literature for more details. The focus 
lies on the description and the testing of the new code elements. 
These are mainly improvements of the hydrodynamics part that 
enhance the accuracy by a careful accounting of the so-called 
"grad-h"-terms and the inclusion of magnetic field evolution. We 
present a formulation of the magne tic field evo l ution in ter ms of 
the so-called Euler potentials (e.g. lEuledll769l: ISternlll994l) that 
are advected with the flow and thus guarantee the MHD flux- 
freezing condition. For comparison, our code also allows to evolve 
the B-fields via a SPH-di scretized version of the MHD-equations 
jPrice & Monaghanir2005l) . 

The paper is organized as follows. Sectionf2]describes the details of 
the numerical methods and their implementation. It includes a brief 
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summary of code elements described previously (Sec. l2.lt . a con- 
cise summary of the hydrodynamics plus gravity as derived from a 
Lagrangian accounting for the so-called "grad-h"-terms (Sec. 12.2b . 
and a detailed description of the treatment of the magnetic field 
(Sec. l2.3t . Various tests of the different code elements are presented 
Sec. [3] The paper is summarized in Sec. [4] 



parallelized and allows in its current form the simulation of sev- 
eral million particles on a medium-sized (24 processor) shared- 
memory computer. Forces emerging from the emission of grav- 
itational waves are treated in a simple a ppro ximation. For more 
details , we refer to lRosswog et alj bOOOl) and iRosswog & Daviej 
d2002h . 



2 METHOD DESCRIPTION 

This section describes the methods and the implementation 
of the various elements of our new (magneto-)hydrodynamics 
code. Historically, numerical MHD schemes have first been 
developed for grid-based methods. Several implementations of 
mag netic fields into smoothed p article hydrodynamic s exis t, 
e.g. iGingold & Monaghanl dl977h : IPhillips & Monaghanl dl985h : 
lB0rve et alj J200lh . but their success has somewhat been hampered 
by several numerical difficulties, not least of which is the difficulty 
in fulfilling the V ■ B = O-constraint. In grid-based methods vari- 
ous techniques exist to enforce this constraint. Our main reasons for 
choosing the smoothed particle hydrodynamics method are the ex- 
act conservation of mass, energy, linear and angular momentum by 
construction and its ease in treating vacuum. On top of that, a La- 
grangian scheme is a natural choice to simulate the highly variable 
geometry that occurs during stellar collisions. The SPH meth od has 
been d e scribed and review e d many times in th e literature (e.g. lBenzl 
( 1990): lMonaghanl dl992h : lMonaghad J2005h ) and this will not be 
repeated here. In comparison with our earlier work we have made 
substantial changes in the implementation of the hydrodynamics 
equations that are documented and tested below, see Sec. 12.21 and 
13.21 In summary, although the new and more sophisticated formu- 
lations improve the accuracy slightly, the changes in the results are 
only minor for the applications that we discuss here. 
The other major change is the inclusion of magnetic fields, which 
is described in Sec.|2.3|and tested in Sec. 13. 31 



2.1 Summary of code ingredients described elsewhere 

The neutron star matter is modeled with th e temperatu r e-depen dent 
relativistic mean-field equation of state of IShen et alj d 1998al fbh. It 
can handle temperatures from to 100 MeV, electron fractions 
from Y e = up to 0.56 and densities from about 10 to more than 
10 15 g cm -3 . No attempt is made to include matter constituents 
that are more exotic than neutrons and protons at high densities. For 
more details on this topic we refer to lRosswog & Daviesl d2002l) . 
The code also contains a detailed multi-flavor neutrino leakage 
scheme. An additional mesh is used to calculate the neutrino opac- 
ities that are needed for the neutrino emission rates at each particle 
position. The free-streaming and neutrino diffusion limit are repro- 
duced correctly, the semi-transparent regime is treated by an inter- 
polation between these limiting cases. The neutrino emission rates 
calculated in this way are used to account for the local cooling and 
the compositional changes due to weak interactions such as elec- 
tron captur es. A detailed description of the neutrino treatment can 
be found in lRosswog & L iebendorferl ( 120031) . 
At present, the self-gravity of the fluid is treated in a Newtonian 
fashion. Both the gravitational forces and the search for the par- 
ticle neighbors that are required, for example, to calculate densi- 
ties or pressure gradients, are performed with a binary tree that is 
based on the one described in iBenz et alj H990I) . These tasks are 
the computationally most expensive part of the simulations and in 
practice they completely dominate the CPU-time usage. The tree is 



2.2 Hydrodynamics 

We are interested in a numerical solution of the Lagrangian, self- 
gravitating Euler equations of ideal hydrodynamics: 

at p 

where if is the fluid velocity, p the pressure, p the mass density and 
$ the gravitational potential. The evolution equation for the specific 
internal energy, u, is 

— = V • v (2) 

at p 

and the density evolves according to 
dp 



dt 



-pV • v. 



(3) 



2.2.1 SPH with "grad-h" -terms 



The exact conservation of energy, linear and angular momentum 
even in the discretized form of the equations is a major strength 
of smoothed particle hydrodynamics. The equations of motion can 
be derived by using nothing more than a suitable Lagrangian, the 
first law of thermodynamics and a prescription on how to obtain an 
density estimate via summation. 

The first derivation of the SPH-equations that takes terms from 
the derivatives of the smooth ing lengths into account goes back to 
Nelson & Papaloizou ( 1994). More recently, Springel & Hernquist 
j2002h and [Monaghanl J2002T) derived the corresponding equations 
from a Lagrangian in two different ways. The equations of this 
more recent approach are less cumbersome to implement, but they 
need an extra iteration for each particle at each time step to make 
the density and smoothing length consistent with each other. The 
advantage of a derivation from a Lagrangian is -apart from its 
elegance- that the resulting equations guarantee the conservation of 
the physically conserved quantities, provided that the Lagrangian 
possesses the right symmetry propert ies. An in-d epth analysis of 
various SPH-variants can be found in |PriceU2004l) . Without going 
into details of the derivations, we will briefly sketch how to arrive 
at the SPH-equations including the so-called grad-h terms. 
The Lagrangian of a perfect fluid is given by jEckartll960T) 



^hyd — 



y - u{p,s] 



dV, 



(4) 



where p is the mass density, v the fluid velocity, u the specific 
energy ("energy per mass") and s the specific entropy. In SPH- 
discretization the Lagrangian reads 
..2 



isPH.h = J]] rrij 



(5) 



where the indexed quantities refer to the values at the SPH particle 
positions. This Lagrangian does not include gravity, the gravita- 
tional terms will be discussed in Sec |2.2.2l The discretized momen- 
tum equation is then found by applying the Euler-Lagrange equa- 
tions 
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d ( 9ispH,h 



-( 

dt \ 



dili 



dL s 



dxi 



= 0, 



(6) 



where Xi and Vi refer to the position and velocity of particle i. The 
first term in Eq. ^ provides the change of the particle momentum, 
rriiVi, the second term in the Lagrangian acts like a potential. The 
second term in Eq. l[6ll becomes 



dxi 



duj(pj,Sj) 



dxi 



3 



dpi 



P-(7) 

dxi 



The derivative with respect to density can be expressed via the first 
law of thermodynamics, which reads for the adiabatic case 



P 

du = -it dp, 

p 2 

where P is the gas pressure. Therefore, 



du 3 
dp 3 



1 j 



and the momentum equation becomes 



dvi 



3 



. Pj 9 Pj 
' p] dxi 



(8) 



(9) 



(10) 



Eq. {8]l also provides us with the evolution equation for the specific 
energy 

dm _ Pi dpi 
~~dt ~ pi dt ' 



(11) 



For the explicit SPH equations we need to specify a prescription for 



the density and to calculate its derivatives ^s 2 - 
and i ll It . For the density we use 

Pi = mjW(nj,hi). 



and 



dpi 



, see 



Eqs. 



(12) 



Here, W is the SPH smoothing kernel, tv, = \f\j\ — \f\ — rj\ 
and hi is the smoothing length of particle i. Throughout this paper 
we use the standard cubic s pline kernel most often used in SPH 
(Monaghan & Lattanzicl l 19851) . Note that contrary to some earlier 
formulations of SPH, only the smoothing length of the particle it- 
self, hi, is used rather than some average. To obtain adaptivity, we 
determine the smoothing length evolution from the density (which 
for equal mass particles is equivalent to a dependence on the parti- 
cle number density). In 3D we use 



, , mi 

hi=ri\ — 

Pi 



1/3 



(13) 



where r\ is a parameter typically in a range between 1.2 and 1.5. A 
careful discussion of the choice of r\ can be found in Price (2004). 
The density pi depends on hi, see Eq. d 1 2b . and vice versa, see 
Eq. d 1 3h . so an iteration is required to reach consistency. Typically 
we iterate until the relative change between two iterations is smaller 
than 10~ 3 . 

By straight forward differentiation of the density sum, Eq. dl2t . one 
obtains 



~dT = f[:^Z m i^ ViWi oi h i 



where Vi 

dPj = _ 

dxi S 



EdWjk(hj) 

k 



(14) 



(15) 



where 



n k = [i-^L.y m " Wkl {h k ) 

\ dp k dh k 



(16) 



With the derivatives, Eqs. J14t and 1151 . at hand the energy equa- 
tion, Eq. dl lb . becomes 



du i]h 1 Pi yr-^ 

/ 

dt Qi p\ ^ 



ijVijViWij(hi 



(17) 



and the momentum equation reads 

dvj^ _ \ - f Pi 
dt ~~ ~ ^ 



m i \ TruViWijihi) + —tV.U;, (./',! } -(18) 
\ilpi "jpj 



We calculate the density via summation, see Eq. l !12t . which solves 
the continuity equation without the need to explicitely evolve the 
density. 

For reasons of reference we provide the SPH equations with aver- 
aged smoothing lengths, h%j = (hi + hj ) /2, that are still widely 
used and that we have used in earlier calculations (this will be re- 
ferred to as the "old equations" or the htj -version). In this /in- 
version the density summation reads 



(Pi)ij = 2_^ m 3 W (\ri-r J \,hi :j )\), 



the energy equation is 
(~^) ■ ■ = m j v ij ViW tJ (h ij ) 

3 ' 3 

and the momentum equation reads 



(19) 



(20) 



=-E^{S + S K " ,7 ' 1 i2l] 

l 3 j 



2.2.2 Self- gravity and gravitational softening 

Most often gravitational softening is done by -physically 
motivated- but still ad hoc recipes. It is, however, possible to de- 
rive the gravitational softening terms self-consistently from a La- 
grangian and to also take the effects from a locally varying smooth- 
ing length into account jPrice & Mon aghan 2007), similar to the 
case of the hydrodynamics equations. 

If gravity is taken into account, a gravitational part has to be added 
to the Lagrangian, Lsph = £sph,1i + ^SPH.g- This gravitational 
part of the Lagrangian reads 



ispH, g = - J~] mj$j 



(22) 



where <3?j is the potential at the particle position j, 3>(rj). The po- 
tential $ can be written as a sum over particle contributions 



$(f) = -G^m J 0(|f-r J |,/ l ), 



(23) 



where h is the smoothing length, <j> the gravitational softening ker- 
nel, and the potential is related to the matter density by Poisson's 
equation 



V 2 $ = 4nGp. 



(24) 
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If we insert the sum representations of both the potential, Eq. {23}, 
and the density, Eq. into the Poisson equation, Eq. d24b . we 
obtain a relationship between the gravitational softening kernel, <f>, 
and the SPH-smoothing kernel W: 

W(\r>-? i \,h) = -±JL (r 2 ^(\r-r 3 \,h)) . (25) 

Here we have used that both <j> and W depend only radially on the 
position coordinate. 

Applying the Euler-Lagrange equations, Eq. (HQl , to Lsph k yield s 
the particle acceleration due to gravity ( Price & Monagh arjj |2007h 



dvi 



dt 



- G £< 

3 

2 ^ 



4>i 3 {hi) + ck'ijQij) 



,(26) 



where ^ = d(f)/d\fi ~ fj\ and iij = i^* ■ The first term in 
Eq. l !26t is the gravitational force term usually used in SPH. The 
second term is due to gradients in the smoothing lengths and con- 
tains the quantities 



dh k \ ^ ^ d(j>kj{h k ) 



dh k 



(27) 



and the defined in Eq. <16t . Formally, it looks very similar to 
the pressure gradient terms in Eq. J 1 8b with correspond- 
ing to Pfc / ' p\. As ffc is a negative definite quantity, these adaptive 
softening terms act against the gas pressure and therefore tend to 
increase the gravitational forces. 

The explicit forms of (f), <f>' and d<j>/ dh for our cubic spl i ne ker - 
nel can be found in Appendix A of iPrice & Monaghanl [2007). 
The gravitational softening procedure obviously only applies to 
int eracting SPH p articles. Generally, we use a binary tree based 
on iBenz eta l. (1990) for the long-range part of the gravitational 
forces. Depending on the choice of the tree opening parameter, 9, 
for each particle a list of nodes is returned whose gravitational in- 
fluences are calculated up to quadrupole order. Forces from nearby, 
interacting particles are calculated via direct summation according 
to the above prescription. 



2.2.3 Dissipation 

The conservation of mass, energy and momentum across a shock 
front requires kinetic energy to be transformed into internal energy. 
Physically, this transformation is mediated via viscosity and usu- 
ally occurs over a length scale of a few mean free paths in the gas. 
This length scale is generally much shorter than any numerically 
resolvable length and thus in the numerical discretization the 
transition appears to be a discontinuity. There are two approaches 
to treat this problem, either by solving a local Riemann-problem 
as in Godunov-type methods, or by adding a controlled amount 
of viscosity artificially to broaden the shock to a numerically 
resolvable width. Not doing so results in unphysical oscillations in 
the post-shock region. While the first approach is certainly more 
elegant, the second one is more robust and offers advantages in 
cases where the analytical solution to the Riemann problem is not 
known, for example, in the case of a complicated equation of state. 
Usually, the artificial viscosity approach is used in SPH and shock 
fronts are usually spread across a few smoothing lengths (rather 
than a few mean free paths) to make them numerically treatable. 



We use an artificial viscosity p r escript io n that is oriented 

at Riemann-solvers (Monag hanl 1 19971 : IPrice & Mon aghan 

l2005h together with tim e dependent v i scosit y parameters 
l Morris & MonagharJ 1 19971 ; iRosswog et al] l2000h so that the 
dissipative terms are only applied if they are really necessary to 
resolve a shock. The additional term in the momentum equation 
reads 



dvi 



dt 



\t)Vs 



P'J 



(28) 



where ptj = (pi + Pj)/2 and the symmetrized kernel gradient is 
given by 



1 



±-ViWii(hi) + i-ViWii(ftj) 



(29) 



The signal speed, v s i s , is the fastest velocity with which informa- 
tion can propagate between particle i and j and for the hydrody- 
namic case it is given by 



c t + Cj 



(30) 



where is the sound velocity of particle k. The time dependent 
parameter that controls the amount of artificial dissipation, afj , is 

4 V = \ (M<) + ■ (/< + fi)} > (3D 

where the fk are the switches suggested by Balsara| d 19951) to sup- 
press effects from artificial viscosity in pure shear flows 

IV -H\ h 



fk = 



(32) 



|V • v\ k + |V x v\ k + 10- 4 ■ c k p k /h k ' 

Here the small additive term in the denominator has been inserted 
to avoid the switch from diverging in case both V ■ v\ and | V x v\ 
tend to zero. The dissipative term in the evolution equation of the 
specific energy reads 

mj — < — [vi ■ ey - Vj ■ e». 
pij I 



' dui 



dt 



\ViWii 



It is straight forward to check that the total energy ■& . \m,iv\ + 
rriiUi) — 0, i.e. that the applied dissipative terms are consistent 
with each other and conserve the total energy. 
The viscosity coefficients , en, are calculated accordin g to an addi- 
tional evolution equation (M orris &Monaghanll 19971) 



den 
~~dT 



ao 



+ Si 



where the decay constant is 

hi 

T~i = 

0.2 c, 

and the source term dRosswog et alj200ol) 

Si = max(-V ■ v, 0)(2 - oh) 



(34) 



(35) 



(36) 



is used. In the absence of compression (V • v > 0) the parameter 
on decays to a minimum value ao, which we choose as 0.1. Note 
that this is more than an order of magnitude below the old SPH- 
prescriptions that used values of a few, and it is further reduced by 
the switch, Eq. d32t . In case of compression, V ■ v < 0, a can rise 
to values of up to 2 in the case of strong shocks. 
Under certain circumstances it is desirable to add a small amount 
of thermal conductivity. This leads to an extra term in the evolution 
equation of the specific energy 
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' du. 



dt 



Vsi g a.ij(ui - Uj) 



Pi} 



(37) 



where — (af + af)/2. The conductivity coefficient a k is 
evolved according to 



dt r fc 



(38) 



where the d ecay constant is the same as above and the source term 
is given by dPrice & Monag han 2005) 



=0.1h k \V 2 u k \ 



(39) 



where we use the Brookshaw-type (Brookshaw 1985) second 
derivative 



(V 2 u) =2Y mj Ui - UjlViWij] 



(40) 



2.3 Magnetic field 

The continuum equations to be solved are those of ideal magneto- 
hydrodynamics. The corresponding momentum equation reads 



dv i _ 1 dS ij 
dt p dx 3 ' 

where the stress tensor is given by 



S'- 



-PS 13 + 



Mo V 



B l B 3 



-B'S 



(41) 



(42) 



and the B k are the components of the magnetic field strength. This 
form accounts for B(X7-B) terms which are needed for momentum 
conservation in shocks but on the other hand are the cause of all the 



numerical instability, see lPrice & Monaghan ( 2004a) for a detailed 
discussion. Both the energy and the continuity equation have the 
same form as in pure hydrodynamics, compare Eqs. ((2J and {3). 
Here we present a discretized SPH-formulation including Euler po- 
tentials. This is our method of choice to evolve the magnetic field. 
For comparison and since some of the equations will be needed 
later, we also summarize a more straight forward SPH disc retiza- 
tion of the MHD equations due to lPrice & Monaghan! (2005). Both 
methods are implemented in the code and it is straightforward to 
switch between the two. 



where flj is the variable smoothing length term defined in Eq. (16}. 
The SPH formulation of the Lor entz force follows natura lly from 
the Lagrangian and is given by jPrice & Mon aghan 2004^) 



dv 



i,mag, 1 

~~dt 



e ~ ' ' 

i 



po 



Bi [Bi-ViWijihi)] 



B 3 [Sj-ViWiAhj)} 



,(44) 



where the terms correspond to the isotropic (due to gradients in 
magnetic pressure) and anisotropic magnetic force (due to field line 
tension) respectively. This exactly momentum-conserving form of 
the anisotropic SPMHD force is known to be unstable to a particle 
clumping instability in the regime where the magn etic field is dom - 
inant over the gas pressure (i.e. for tension forces) jMorrisll996b ll 
lB0rve et al]2004l) . Whilst typical magnetic field strengths found in 
compact objects mean that most simulations we will perform will 
lie in the regime where the above formulation is stable, a simple 
solution in the unstable regime is to replace the anisotropic com- 
ponent of the m agnetic force with a formulation that vanishes for 
constant stress (Morris 1996b), given by 



dvi 



mag, 2 



dt 



/.to 



-E 

3 

Em, 
— 
I In 



Bf/2 
^P 2 



B 2 /2 



Bi(Bi ■ X7 t W lj ) 



B ji B 3 



ViWi, 



PiPj 



Using this force means that momentum is no longer conserved ex- 
actly on the anisotropic term, however the effect of this small non- 
conservation on shocks (w here good conservation is critical) proves 
minimal (see lPricd J2004h ). 



2.3.2 Dissipation 

Dissipation te rms necessary for the treatm ent of MHD shocks were 
formulated by Price & Monaghan] d2004al) . The induction equation 
contains a dissipative term corresponding to an artificial resistiv- 
ity, ensuring that strong gradients in the magnetic field (i.e. current 
sheets) are resolved by the code. This term is given by 



d_ 

dt 



E 



Pij 



Bi)\VWi; 



(46) 



(45) 



2.3.1 Smoothed Particle Magnetohydrodynamics 

The following SPH discretization goes back to 
( Phillips & Monaghan 1985), a nd has been e xtended and refined 
recent ly by iPrice & Monaghan! d2004a l |(3) and IPrice & M onaghan 
(2005). This algorithm has been extensively tested on a wide range 
of problems used to benchmark grid-based MHD codes. As in the 
hydrodynamic case, the formulation can b e elegantly derived from 
a Lagrangian JPrice & Monaghanl l2004bl) . guaranteeing the exact 
conservation of energy, entropy and momentum. The "grad-/i" 
for mulation of the SPM HD equations was derived in this manner 
by|p rice & Monaghanl |2004b) and we use this formulation here. 
The magnetic flux per unit mass B/p evolves according to 



d 
di 



. Pi . 



P 2 ^ 



E 



TTljVijBi 



ViW tj {hi), 



(43) 



where the energy equation contains a corresponding term 

af j (t)v stg 



\dt) dis ,B ^ 



Mop? 



-(Bi - BjYIVWi: 



(47) 



It is straightforward to demonstrate that this term gives a pos itive 
definite contribution to the entropy dPrice & Mon aghan 2004ah. 

In the magnetic field case we use a simple generalization of 
the signal velocity, Eq. J30b . given by 

Vsig = i(Vi + Vj) - {Vij ■ Sij), (48) 

where v is the maximum propagation speed for MHD waves given 
by 

1/2 



1 

71 



fcf + «a) + 



Popi 



,(49) 
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where v\ 



— *- is the Alfven speed. 



2.3.3 Euler potentials 

A key problem associated with the simulation of MHD phenomena 
is the maintenance of the divergence-free condition associated with 
the magnetic field. Whilst various methods for correcting the field 
produced by the standard SPMHD evolution Eq. d43t are possible 
l lPrice & Monaghanl 12005ft . we can avoid the problem entirely by 
formulating the magnetic field such that the divergence constraint 
is satisfied by construction. Use of the magnetic vector potential 
is one such construction. However for particle methods a natural 
choice is the so-c alled 'Euler p otentials' (originally formulated by 
lEulerld 1769b - see Stern ( 1970)) but also referred to as the 'Clebsch 
formulation' (e.g. [Phillips & Monaghaiil l 1985b . In this formulation 
the magnetic field is represented as 



B — Vq x V/3. 



(50) 



Geometrically, the Eu ler potentia ls can be thought of as magnetic 
field line labels (e.g. ISternlll966l) : the magnetic field lines corre- 
spond to the intersections of surfaces of constant a with surfaces of 
constant (3, see Fig.Q] 

The Euler potentials can be easily related to a vector potential 
which can be of the form 



A = aV(3 + V£ 



A = -/3Vq + V^, 



(51) 



(52) 



where £ and tp are arbitrary smooth functions. It is straight forward 
to show that these vector potentials yield the B -field: V x A = 
Va x V/9 = B. As the Euler potentials only contain two indepen- 
dent variables (rather than the three components of A), they corre- 
spond to an implicit choice of a gauge for the vector potential and 
that is maintained exactly during the further evolution. Taking the 
divergence of Eq. J50b demonstrates that the V • B = constraint 
is satisfied by construction. 

The condition that the magnetic field is frozen in translates into a 
pure advection of the Euler potentials with the particles: 



= 0. 



(53) 



den _ dfh 
dt ~ ' dt 

This advection property Eq. J531 > means that the evolution of an ar- 
bitrary magnetic field can, in principle, be reconstructed from a hy- 
drodynamic simulation given the initial and final particle positions 
(as long as the feedback from magnetic forces does not change the 
flow). 

There are, however, some restrictions of the Euler potentials in 
comparison to a MHD scheme where all three components of the 
magnetic field are evolved. These are: 

i) the calculation of the force involves second derivatives of the 
potentials, which may be less accurate. 

ii) zerc0 dissipation (i.e. no reconnection of field lines) 

iii) non-linearity of initial conditions - that is, for a given B it is 
a non-trivial task to obtain the corresponding Euler potentials. 



1 In our scheme there is some smoothing of the Euler potential gradients 
in the calculation of B from our use of Eqs. {59} and (60} 



P= const. 




oc= const. 



Figure 1. Euler potentials: intersections of surfaces of constant values a 
and j3 correspond to magnetic field lines. 



With regards to i) the tests presented here using the Euler potential 
formalism show no significant differences in force accuracy com- 
pared with similar tests shown using the standard SPMHD formal- 
ism. With regards to ii) we will demonstrate below how this re- 
striction can be overcome by a simple modification to the Euler po- 
tential evolution. With regards to iii) the non-linearity of the Euler 
potentials may present some difficulty for the setting up of compli- 
cated initial conditions, but given the uncertainty of magnetic field 
configurations in most compact objects, we choose a simple initial 
configuration, so it does not present an immediate stumbling block 
for our simulations. 

In two dimensions the Euler potentials are equivalent to a vector 
potential formulation with a = A z and f3 = z. We calculate the 
gradients of the Euler potentials in a way so that gradie nts of linear 
functions are reproduced exactly, see e.g. lPricel 1 20041) . The gradi- 
ent of the product of the density and an arbitrary quantity A reads 
in standard SPH discretization 



[V{pA)]i 



m 3 A 3 ViWij{hi) 



(54) 



If we use the RHS of this equation on both to the left and the right 
side of the equal sign and insert on the right the Taylor expansion 
of Ai around ft 



A, 



Ai + (rj - fif 



OA 

Qfp- 



one finds 



3 



OA 
<9rf 



E 



(55) 



m.i^-rA^Wij, (56) 



where we use Greek letters as summation indices, Latin ones for 
the particle identities and the kernel gradients are evaluated with the 
smoothing length hi. This equation can be solved for the gradient 
of A at the position f t : 



dA 



3 



where the matrix is given by 

3 



(57) 



(58) 
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and V£' is the /^-component of the gradient evaluated at position r,:. 
Applied to the Euler potentials this yields 

(V"a)i = (xTr^msfa -aJVlWiiihi), (59) 

3 

This formulation involves the inversion of a 3 x 3-matrix, x< f° r 
each particle. This can be done analytically and the matrix only 
needs to be stored for one particle at a time. The summations in 
Eqs. ([59} and d60b do not involve densities, therefore they can be 
conveniently be calculated in the density loop for subsequent use 
in the force calculation. 

Whilst in principle it is possible to formulate the magnetic forces 
using direct second derivatives of the Euler potentials - e.g. mak- 



ing us e of the SPH second derivative formul ations of Brookshaw 



1985) generalized to vector derivatives by lEspanol & Revenga 



2003)- it is not possible to do so and at the same time maintain the 



conservation of linear momentum in the force formulation, as the 
force involves a combination of first and second derivatives of the 
potentials which cannot be symmetrized. For this reason we simply 
use the usual force Eq. J44b or Eq. ([45} where B is the magnetic 
field computed using Eqs. ( 1591 ), d60t and ([50}. The tests presented 
in $ 3. 3 1 demonstrate that the resulting force is no less accurate than 
when the induction equation is used to evolve the magn etic field. 
A similar conclusion was reached bv lWatkins eTafl Jl996h . who, in 
formulating Navier-Stokes type viscosity terms for SPH, found that 
using a nested fir s t deri vative could in fact be more accurate than 
using lBrookshawUl985l) type terms. 



2.3.4 Euler potentials with dissipation 

The standard advection of the Euler potentials with the SPH par- 
ticles results in zero dissipation of magnetic field lines and no re- 
connection. However, for problems involving shocks it is necessary 
to add dissipative terms that make the discontinuities numerically 
treatable by spreading them over a few smoothing lengths. In order 
to do so we propose a simple modifica tion of the Euler potential 
evolution based on the Monagha^ 1 1997b formulation of SPH dis- 
sipative terms 



(dm \ 

\ dt J diss 



\ dt J diss 



i 



Pij 
Pij 



(ai-a^lVWyl, (61) 



GSi-ftJIVWyl, (62) 



which are SPH representations of the equations 

d(3 



da 2 
dt 1 



dt 



r?V 2 /3. 



(63) 



where 77 ~ a i: jV S igh. Rigorous conservation of energy would re- 
quire computing the resulting change in dB/dt according to 



dt 



\ dt J \dtj 



(64) 



which should be computed using the SPH formulations Eqs. d59t 
and d60t for the gradients. However this would require an additional 
pass over the particles to compute the terms Eqs. d61t and ( 1621 ) (af- 
ter the density summation) before substituting the result into the 



SPH expression for Eq. d64t . An approximate, but much more ef- 
ficient solution is to simply compute the energy input according to 
Eq. d47| > using the B calculated from the Euler potentials. For the 
shock tube tests we find that this approximate approach is more 
than satisfactory. 

Our sole purpose in formulating dissipative terms for the Eu- 
ler potentials is to provide a mechanism to ensure that discontinu- 
ities in the magnetic field are treated appropriately by the numerical 
method (i.e. resolved over a few smoothing lengths). As such the 
terms given above do not, and are not intended to, correspond to a 
rigorous formulation of Ohmic dissipation using the Euler poten- 
tials. Indeed a more detailed derivation demonstrates that, whilst 
such a formulation would include terms of the form Eq. d63t . ad- 
ditional terms would also be required in order for the dissipation 
to correspond meaningfully to the usual Ohmic dissipation terms 
added to the MHD induction equation. 



2.4 Time integration 

The calculation of the gravitational forces is -together with the 
neighbor search- the computationally most expensive task. It is 
therefore advantageous to choose a time integration method that 
only requires one force evaluation per time step. 
The basic integration scheme that we use is a s econd order accurat e 
MacCormack predictor-corrector method (e.g. lLomax et al]|200ll) . 
The predictor is given by 



2/i+i 



yt + AUy'i, 



(65) 



the corrector is 

Vi+i =Vi + At» (j/j + (66) 

Here, i labels the time step, the primes denote derivatives with re- 
spect to time, y' i+ i denotes the derivatives at the predicted position 
and Ati the used time step. 

Our integration scheme allows for individual time steps, i.e. each 
particle is evolved on its own timestep while the forces at any given 
point in time, t, are calculated from the particle properties interpo- 
lated to t. At the beginning of the simulation, the "desired" time 
step of each particle i is determined 



Ati,des = 0.2 min (At/,i, At c ,i) 



(67) 



where Atf,i — \J hi/\fi\ and fi is the particle's acceleration. 
The Courant-type time step is given by At c ,i = min., ■■ (hj / v B i s ,j) 
where j runs over all neigbors and the particle i itself. At the be- 
ginning of the simulation all particles start out with the same time 
step, dto, that is the maximum time step that fulfills the condition 
dto < mini(Ati,dcs) and the condition that an integer number of 
these time steps is equal to the time of the next data output. In prac- 
tice this means that the particles are running on time steps that are 
slightly smaller than the desired ones, Eq. d67t . During the further 
evolution, we allow the particles to reduce their time step by a fac- 
tor of 2~ n , where 71 is an integer, or, if the new desired time step 
is larger than twice the previously used time step, to increase the 
time step by a factor of 2. The latter is only allowed if the next 
data output time can be hit exactly. Whenever a particle needs to be 
updated, all other particle properties are calculated at this point of 
time by interpolation to obtain the required derivatives. 
The gain in computing time depends to a large extent on the ap- 
plication. In the context of a neutron star merger the gain is only 
moderate, about a factor of two. This is a consequence of the nearly 
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Figure 2. Particle distribution in a slice of 0.5 code units thickness in the 
XY plane of stars with 50000 particles. The left panel shows the stretched 
cubic lattice, the right one the stretched close-packed configuration. 



incompressible neutron star matter that results in a rather flat den- 
sity distribution within the stars. Therefore, the bulk of the particles 
has to be evolved on the shortest time step. In other applications, 
however, say in the tidal disruption of a star by a black hole (e.g. 
iRosswod fcoOS). the gain can be easily more than two orders of 
magnitude. 



3 TESTS 

In this section we will describe tests of the new code elements. We 
start with a description the initial particle setup and discuss why we 
use exclusively equal mass SPH particles for neutron star simula- 
tions. We then present tests of both the hydro- and the magnetohy- 
drodynamics ingredients in one two and three dimensions, where 
the one and two dimensional tests are included as tests of the algo- 
rithms for comparison with other codes and the three dimensional 
tests are performed with the code itself. 



3.1 Particle setup 

It is important to start a simulation from a SPH particle config- 
uration that has been 'relaxed' into its optimal, minimal-energy 
configuration. This is usually done with the full hydrodynamics 
code by applying an additional velocity dependent artificial 
acceleration term fa oc v%. We discuss two different particle setups, 
a cubic lattice and a close-packed configuration, in each case we 
only use particles of the same mass. The case with unequal masses 
will be discussed below. 

The first step is to solve the stellar structure equations in ID 
to find the equilibrium profiles p 1D (r), Y} D (r) and T 1D (r) 
for a neutron star of a specified mass, M ns . Here, p, Y e and T 
are density, electron fraction and temperature. In the next step 
the desired number of particles, N, is distributed inside a unit 
sphere, either on a cubic lattice or as a close-packed configuration. 
To keep the particle mass constant, the number density of the 
particles, n(r), has to reflect the density distribution of the star, 
p 1D (r) — n(r)m, where m = M nB /N is the mass of each SPH 
particle. Subsequently the unit sphere is stretched to the size of the 
neutron star so that the above condition is met. The configuration 
constructed in this way is very close to hydrostatic equilibrium. 
Examples with both types of setups are shown in Fig. [2] To find 
the true numerical equilibrium state we relax this configuration 
with the full hydrodynamics code by apply ing an art ificial, 
velocity dependent damping force (e.g. lRosswog et all 2004'). This 
procedure yields numerical equilibrium conditions with a minimal 
computational effort. To demonstrate that the particle distribution 




Figure 3. Density profile of a 1.4 M© neutron star. The solid red line is the 
result of solving the ID stellar structure equations, the dots are the results 
obtained with the full SPH code at three different resolutions (10 4 , 10 5 and 
10 6 SPH particles). 



really settles to the correct result, we show in Fig. [3] the density 
profile of a 1 .4 Mq neutron star as obtained by solving the stellar 
structure equations in ID ("exact", red solid line). Overlaid are the 
density distributions as obtained by relaxing three neutron stars 
of different numerical resolution: maroon corresponds to 10 000, 
blue to 100 000 and black to 1 000 000 SPH particles. The overall 
agreement with exact result is very good, deviations are only 
visible at the stellar edge where the extreme density gradients are 
challenging. 

As a test of the quality of the initial conditions we set up the 
particles in the described way, but instead of relaxing them, we 
evolve and monitor the kinetic energy that builds up as a result of 
small deviations from the true (numerical) hydrostatic equilibrium. 
We find only very minor differences resulting from the different 
particle setups. For very low particle numbers, say 1 000 particles, 
the gradients in the stellar profile cannot be resolved properly 
and the particles adjust their positions to find the equilibrium. For 
particle numbers in excess of 100 000 the particles smoothly move 
off the initial grid but the overall density structure is practically 
unperturbed. 

Different particle masses are known to introduce numerical noise 
into SPH simulations. While a small range of particle masses 
may be admissible in some applications, we only use equal 
particle masses. As a numerical experiment, we set up a star with 
10 000 particles and a constant particle number density, so that 
the particle masses carry the information of the stellar profile 
p 1D (r) = nm(r). The extreme drop in density towards the 
neutron star surface (caused by the very stiff equation of state) 
translates for this particle setup in particle masses that vary by 
more than a factor of 10 6 . Without further relaxation we let this 
configuration evolve. This worst case setup results in spurious 
particle motions as very light particles are in direct contact with 
heavier particles in the neutron star interior. The slightest noise of 
the heavy particles strongly disturbs the light ones ("ping pong on 
cannon ball effect") and leads to pathological particle densities, 
where low density particles can be found in the interior of the star. 
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Figure 4. Results of the Sod shock tube test in one dimension using 900 
SPH particles setup using unsmoothed initial conditions. Artificial viscosity 
and thermal conductivity are applied to appropriately smooth the shock and 
contact discontinuity respectively. The exact solution is given by the solid 
line. The upper row displays the velocity (left) and the density (right), the 
bottom row shows specific internal energy (left) and the pressure (right). 



Therefore, we only use equal mass particles to keep the numerical 
noise at a minimal level. 



3.2 Hydrodynamics 

3.2.1 ID: Sod's shock tube 

As a standard test of the shock capturing capability we show the re- 
sults of Sod's shock tube test (Sod 1978). To the left of the origin, 
the initial state of the fluid is given by [p, P, v x ]l = [1.0,1.0,0.0] 
whilst to the right of the origin the initial state is [p, P, v x ]r = 
[0.125,0.1,0.0] witti7 = 1.4. The problem is setup using 900 equal 
mass particles in one spatial dimension. Rather than adopting the 
usual practice of s mooth i ng the initial conditions across the discon- 
tinuity, we follow IPricd j2004) in using unsmoothed initial condi- 
tions but applying a small amount of artificial thermal conductiv- 
ity using the switch described in i)2.2.3l The results are shown in 
Fig. [4] where the points represent the SPH particles. For compar- 
ison the exact solution computed using a Riemann solver is given 
by the solid line. 

The shock itself is smoothed by the artificial viscosity term, which 
in this case can be seen to spread the discontinuity over about 5 
smoothing lengths. The contact discontinuity is smoothed by the 
application of artificial thermal conductivity which (in particular) 
eliminates the "wall heating" effect often visible in numerical solu- 
tions to this problem. The exact distribution of particle separations 
in the contact discontinuity seen in Fig. [4] is related to the initial 
particle placement across the discontinuity. 

For this test, applying artificial viscosity and thermal conductivity 
as described, we do not find a large difference between the "grad- 
h" formulation and other variants of SPH based on averages of the 
smoothing length. If anything, the "grad-/i"-terms tend to increase 
the order of the method, which, as in any higher order scheme, 
tends to enhance oscillations which may otherwise be damped, vis- 
ible in Fig.|4]as small "bumps" at the head of the rarefaction wave 
(in the absence of artificial viscosity these bumps appear as small 




Figure 5. Results of the Einfeldt test for the velocity (upper left), density 
(upper right), specific internal energy (lower left) and pressure (lower right). 
400 particles were used in this test. 



but regular oscillations with a wavelength of a few particle spac- 
ings). 

3.2.2 ID: The Einfeldt rarefaction test 

The i nitial conditions of the Einfeldt rarefaction test I Einfe ldt et al.l 
1 19911) do not exhibit any discontinuity in density or pressure, but 
the two halfs of the computational domain move in opposite direc- 
tions and thereby create a region of very low density near the initial 
velocity discontinuity. This low-density region represents a partic- 
ular challenge for some iterative Ri emann solvers which can return 
negative values for pressure/density. Einfel dtet al.l dl99lT) designed 
a series of tests to illustrate this failure mode. The initial condi- 
tions of this test are [p, P, v x ]l = [1.0,0.4,-2.0] for the left state 
and [p, P, v x ]b. = [1.0,0.4,2.0] for the right one. The results from a 
400 particle calculation are shown in Fig.[5]after a time t= 0.2. The 
exact result is reproduced very accurately, only at the fronts of the 
rarefaction waves a small overshoot occurs. The low density region 
does not represent any problem for the method. 

3.2.3 3D: Sedov blast wave test 

In order to demonstrate that our scheme is capable of handling 
strong shocks in three dimensions, we have also tested the code 
on a Sedov blast wave problem both with, see Sec. 13.3.41 and with- 
out magnetic fields. Without magnetic fields the explosion is spher- 
ically symmetric, however for a strong magnetic field the blast 
wave is significantly inhibited perpendicular to the magnetic field 
lines, resulting in a compression along one spatial dimension. Sim- 
ilar tests for both hydrodynamics and MHD have been used by 
many authors - for example by Balsara] d200lh in order to bench- 
m ark an Adaptive Mesh Refin ement (AMR) code for MHD and 
bv lSpringel & Hemquistl fe002h to test a new formulation (entropy 
equation) of SPH. 

The hydrodynamic version is set up as follows: The particles are 
placed in a cubic lattice configuration in a three dimensional do- 
main [—0.5,0.5] x [—0.5,0.5] x [—0.5,0.5] with uniform density 
p — 1 and zero pressure and temperature apart from a small region 
r < R near the origin, where we initialize the pressure using the 
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Figure 6. Results of the hydrodynamic Sedov blast wave test at t = 0.09 at resolutions of 125 000 (top) and 1 million (bottom) particles respectively. The left 
panels show a rendered plot of the density in a slice taken at z = whilst the right panels show the density and radial position of each SPH particle, which 
may be compared to the exact solution given by the solid line. Note that we are showing each individual SPH particle and no averages. 



total blast wave energy E = 1, ie. P = (7 - l)E/(^nR 3 ). We 
set the initial blast radius to the size of a single particle's smooth- 
ing sphere R = 2r/Ax (where 2 is the kernel radius, r;(= 1.5) is 
the smoothing length in units of the average particle spacing as in 
Eq. U3\ and At is the initial particle spacing on the cubic lattice) 
such that the explosion is as close to point-like as resolution allows. 
Boundaries are not important for this problem, however we use pe- 
riodic boundary conditions to ensure that the particle distribution 
remains smooth at the edges of the domain. 

The results are shown in Fig. [6] at t = 0.09. We have used a res- 
olution of 50 3 and 100 3 particles (ie. 125 000 and 1 million parti- 
cles respectively) and we have plotted (left panels) the density in a 
z — cross section slice and (right panels) the density and radial 
position of each particle (dots) together with the exact self-similar 
Sedov solution (solid line). 

We found that the key to an accurate simulation of this problem in 
SPH is to incorporate an artificial thermal conductivity term due to 
the huge initial discontinuity in thermal energy. The importance of 
su ch a te rm for shock problems in SPH has been discussed recently 
bv lPricel ^20041) . In the absence of this term the particle distribution 
quickly becomes disordered around the shock front and the radial 
profile appears to be noisy. From Fig.|6]we see that at a resolution 
of 1 million particles the highest density in the shock at t = 0.09 
is pmax = 2.67 whereas for the lower resolution run p max = 2.1, 
consistent with a factor of 2 change in smoothing length. Using 
this we can estimate that a resolution of ~ 345 3 = 41 million par- 
ticles is required to fully resolve the density jump in this problem 
in three dimensions. Note that the minimum density obtained in the 
post-shock rarefaction also decreases with resolution. Some small- 
amplitude post-shock oscillations are visible in the solution which 



we attribute to interaction of the spherical blast wave with particles 
in the surrounding medium initially placed on a regular (Cartesian) 
cubic lattice. 



3.2.4 3D: Radial oscillation of a neutron star 

As a further test case, we consider the ra dial oscillations of a neu- 
tron star using the Shen equation of state dShenetal.ll9983lbh . The 
initial conditions are a neutron star relaxed into hydrostatic equilib- 
rium as described previously ( i]3.1K given an initial perturbation in 
velocity of the form v = nor where Do is an arbitrary but small am- 
plitude (we choose vo = 0.01c). No artificial viscosity or damping 
is applied for this problem since no shocks are involved. We com- 
pute the problem at low resolution using only 10 000 particles in 
the neutron star. 

The results of this test are given in Fig. [7] which shows the results of 
an integration for 10 oscillation periods, where top and bottom pan- 
els show the total and gravitational potential energy respectively. 
From this Figure it may be observed that the amplitude is main- 
tained almost exactly by the code over the 10 oscillation periods 
(P ~ 0.33 ms) simulated. The residual fluctuations in total energy 
are directly attributable to a combination of the tree opening cri- 
terion (here 9 = 0.9), which we find controls the level of "noise" 
in the total energy curve, and the timestepping accuracy (Courant 
number, here C CO ur = 0.2), which affects the mean curve. 
The fact that our SPH code, even at low resolution, is capable of 
following the neutron star oscillations for many periods without 
significant damping suggests that the code may be an ideal tool 
for studying neutron star oscillation modes. A similar study has 
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Figure 7. Total (top) and gravitational potential energy (bottom) for 10 ra- 
dial oscillations of an initially hydrostatic neutron star using the Shen equa- 
tion of state. 
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Figure 9. Evolution of total energy and angular momentum (both in code 
units) during the orbital revolution of a binary system with twice 1 .4 Mq . 
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Figure 8. Evolution of a neutron star binary system set on a circular orbit 
with a initial separation of ao =52.5 km. Displayed are the orbital separa- 
tion (upper line) and the maximum matter density (lower line). The binary 
system is evolved for as long as 63 orbital periods or about 920 neutron star 
dynamical time scales. 



recently been performed by Monaghan & Price ( 2006) who com- 
pared SPH s imulations of the oscillatio n modes of two dimensional 
"Toy stars" (Monaghan & Price 2004) with exact and perturbation 
solutions, finding good agreement between the two. 



3.2.5 3D: Binary orbit 

As a further test in 3D we set up a neutron star binary system on 
a stable circular orbit and follow its long-term evolution. As initial 
separation we choose a ( ) = 52.5 km (= 35 code units), no gravita- 
tional backreaction forces are applied. To demonstrate that even at a 
very low resolution stable and accurate orbital evolution can be ob- 
tained, we model each neutron star with 1 300 SPH particles only. 
We relax two neutron stars in a corotating frame as described in 



iRosswog et al.l |2004). After a tidally locked equilibrium configu- 
ration has been reached, the velocities are transformed to the space- 
fixed frame and the orbital evolution is followed with the full code 
for as long as 63 orbital periods or approximately 920 dynamical 
time scales of the neutron stars. The evolution of the orbital sepa- 
ration of both neutron stars together with the maximum density in 
the binary system are shown in Fig. [8] The binary stays nearly per- 
fectly on the intended orbit. Due to the finite relaxation time very 
small scale oscillations occur which lead to an exchange between 
orbital and oscillation energy of the stars. This leads to small os- 
cillations of the orbital separation around the initial value. But the 
corresponding deviations are very small (5a max /ao ~ 0.005) and 
they do not grow during the very long evolution time. The central 
density is free of any visible oscillation. 

The evolution of the corresponding total energy, E to t, and the total 
angular momentum, L to t, are shown in Fig. [9] Both quantities are 
excellently conserved: 6E to t/ 'E to t,o < 10~ 3 and SL to t/L to t,o < 
10" 4 . 



3.2.6 3D: Stellar head-on collision 



It has long been known dHernquistll 19931) that using the SPH equa- 
tions derived under the assumption of constant smoothing lengths, 
e.g. in the conventional hij .-formulation summarized in Sec. 12. 2. 1~| 
but still allowing the smoothing lengths to change in practice, can 
in extreme cases lead to substantial er rors in the conservation of 
energy. For example, iHemquistl ( fl993l) found a non-conservation 
of energy on a ~ 10-%-level for a violent head-on collision of two 
polytropic stars. To quantify this non-conservation for MAGMA 
we perform a similar head-on collision between two neutron stars. 
Two neutron stars (35 000 particles each) of 1 .4 Mq obeying the 
IShenetal.l ll998b) equation of state are set to an initial separation 
of 35 code units (=52.5 km) and provided with an initial relative 
velocity of 0.1 c. Fig. 1 101 shows two density snapshots with over- 
laid velocity field during the evolution. 

In FigQT]we show the evolution of the total energy both for the hij- 
formulation and the new "grad-/T,"-versi on (both are normalized to 
their initial values). As in previous work dRosswog & Daviesl2002t 
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Figure 10. Head-on collision of two neutron stars. Such a head-on collision of two stars is considered a worst-case situation for the non-conservation of energy 
due to changes in the smoothing lengths, see text for details. 
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Figure 11. Comparison of the energy conservation of both SPH- 
formulations (see text for details) for the above head-on collision. Both 
formulations yield very good conservation properties. The dominant error 
sources are the tree opening criterion and the time stepping criterion. To 
see a difference from the grad-h-terms, the tree opening criterion had to be 
reduced to = 0.1 and the pre-factor in Eq. |67| to 0.1. 



iRosswog & Liebendor feil r2003l ; iRosswog et alj [2003) we use on 
average 100 neighbors for the hij -formulation, for the "grad-/i"- 
version we use a constant of tj = 1.5 in Eq. d 1 3 b - 
Generally the energy is very well conserved in both cases and the 
non-conservation is determined by the tree-opening criterion and 
the time stepping accuracy. To see a difference between both for- 
mulations, we reduced the tree opening criterion to 6 — 0.1 and 
the pre-factor in Eq. l |67t to 0. 1 . Both codes conserve energy in this 
challenging problem to better than about 10 -3 with the "grad-/i"- 
version showing a slightly better performance. As in the other tests 



presented here, wee see a small improvement, but no major change 
due to the use of the grad-h,-terms. 



3.3 Magnetohydrodynamics 

3.3.1 ID: Brio-Wu shock tube test 



The magnetic shock tube test of iBrio & Wi] dl988l) has be- 
come a standard test case for numerical MHD schemes that has 
been widely used by many authors to benchmark (mainly grid- 



based) MHD codes (e.g.lStone et al.ll992: Dai & Woodwar d! 1994 



iRvu&Jonesll 19951 : Balsar; 3 119981) . The Brio-Wu shock test is the 
MHD analogon to Sod's shock tube problem that was described 
in Sec. 13.2.11 but here no analytical solution is known. The MHD 
Riemann problem allows for much more complex solutions than 
the hydrodynamic case which can occur because of the three dif- 
ferent types of waves (i.e. slow, fast and Alfven, compared to just 
the sound waves in hydrodynamics). In the Brio-Wu shock test the 
solution contains the following components (from left to right in 
Fig.ll2t: a fast rarefaction fan and a slow compound wave consist- 
ing of a slow rarefaction attached to a slow shock (moving to the 
left) and a contact discontinuity, a slow shock and a fast rarefaction 
fan (moving to the right). It has been pointed out, however, that 
the stability of the unusual compound wave may be an artifact of 
the restriction of the symmetry to one spatial dimensi on whilst al- 
lowin g the magnetic field to vary in two dimensions dBarmin et al.l 

mil). 

Here we present the first results using the Euler potential formu- 
lation, see o|2.3,4| Results of this problem using Smoothed Par- 
ticle Magnetohydrodynami cs (S PMHD ) have been presented by 
IPrice & MonagharJ j2004al) and IPricd J2004I) . The Euler poten- 
tials show a distinct improvement over the standard SPMHD re- 
sults. The initial conditions on the left side of the discontinuity 
are [p,P,v x ,v y ,B y ]- L = [1, 1, 0, 0, 1] and [p, P, v x , v y , B y ] n = 
[0.125, 0.1, 0, 0, —1] on the right side. The component of the 
magnetic field is B x = 0.75 everywhere and a polytropic exponent 
of 7 = 2.0 is used. Using the Euler potentials the components are 
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Figure 12. Results of the Brio & Wu MHD shock tube test at t = 0.1 using 631 particles and the Euler potential formulation. For comparison the numerical 
solution taken from Balsara 1 1998) is given by the solid line. The solution illustrates the complex shock structures which can be formed due to the different 
wave types in MHD, including in this case a compound wave consisting of a slow shock attached to a rarefaction wave. 



given by a = —B y x (equivalent to the vector potential A z ) and 
P = z (or more specifically V/3 = z) and the B x component is 
treated as an external field which requires adding a source term to 
the evolution equation for a as discussed in ^2.3.41 Particles are re- 
stricted to move in one spatial dimension only, whilst the magnetic 
field is allowed to vary in two dimensions (that is, we compute a v y 
but do not use it to move the particles). This is sometimes referred 
to as a "1.5D" approximation. 

We setup the problem using 631 equal mass particles in the domain 
x 6 [—0.5, 0.5] using, as in the hydrodynamic case, purely discon- 
tinuous initial conditions. Artificial viscosity, thermal conductivity 
and resistivity are applied as described in ^2.2.31 and i]2.3.4| The 
results are shown at t = 0.1 in F ig.[72l For comparison the numeri- 
cal solution from Balsara (1998) is given by the solid line (no exact 
solution exists for this problem). The solution is generally well cap- 
tured by our numerical scheme. Two small defects are worth noting. 
The first is that a small offset is visible in the thermal energy - this 
is a result of the small non-conservation introduced by use of the 
Morris formulation of the magnetic force (required for stability, see 
Eq. (145b V Secondly, the rightmost discontinuity is somewhat over- 
smoothed by the artificial resistivity term. We attribute this to the 
fact that the dissipative terms involve simply the maximum signal 
velocity v S i 9 (that is the maximum of all the wave types). Ideally 
each discontinuity should be smoothed taking account of its indi- 
vidual characteristic and corresponding v s i g (as would occur in a 



Godunov-MHD scheme). Increasing the total number of particles 
also decreases the smoothing applied to this wave. 



3.3.2 2D: Current loop advection problem 

A simple test problem for MHD is to compute the advec- 
tion of a weak magne tic field loop. This test, introduced by 
Gardin er & Stond J2005h in the development of the Athena MHD 
codfl presents a challenging problem for grid-based MHD 
schemes requiring careful formulation of the advection terms in the 
MHD equations. For our Lagrangian scheme, this test is straight- 
forward to solve which strongly highlights the advantage of using 
a particle method for MHD in problems where there is significant 

motion with respect to a fixed reference frame. 

We setup the problem here following iGardiner & Stonel d2005l) : 
The computational domain is two dimensional with x £ [—1,1], 
y 6 [—0.5, 0.5] using periodic boundary conditions. Density and 
pressure are uniform with p = 1 and P = 1. The particles are 
laid down in a cubic lattice configuration with velocity initial- 
ized according to v = (vo cos 6, vo sin 0) with cose = 2/VB", 
sin 8 = l/\fE and Vo = 1 such that by t — 1 the field loop will 
have been advected around the computational domain once. The 



http://www.astro.princeton.edu/~jstone/athena.html 
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Figure 13. Magnetic field lines in the current loop advection test, plotted 
at t = (top) and after 1 000 crossings of the computational domain (bot- 
tom). 



magnetic field is two dimensional, initialized using a vector poten- 
tial given by 



A z 



A (R 




r < R, 
r>R, 



(68) 



where Ao = 10~ 3 , R = 0.3 and r — *J x 2 + y 2 . The ratio of 
thermal to magnetic pressure is thus given by /3 p i as = P/(^B 2 ) = 
2 x 10 6 (for r < R) such that t he magnetic field is passively ad- 
verted. iGardiner & Stonel {2005) show the results of this problem 
after two crossings of the computational domain, by which time 
the loop has either been significantly diffused or has disintegrated 
into oscillations depending on details of their particular choice of 
scheme. The advantage of a Lagrangian scheme is that advection is 
computed exactly, and using our Euler potential formulation (which 
in two dimensions is equivalent to a vector potential formulation 
with a — A z and j3 = z) for the magnetic field, this is also true 
for the evolution of the magnetic field. The result is that the field 
loop is adverted without change by our code for as long as one may 
care to compute it. This is demonstrated in Fig.[T3]which shows the 
magnetic field lines at t — (top) and after 1 000 (this is not a mis- 
print!) crossings of the computational domain (bottom), in which 
the field configuration can be seen to be identical to the top fig- 
ure. The magnetic energy ( not sh own) is also maintained exactly, 
whereas Gardine r & Stonel |2005) find of order a 10% reduction in 
magnetic energy after two crossings of the domain. 
In a realistic simulation involving MHD shocks there will be some 
diffusion of the magnetic field introduced by the addition of ar- 
tificial diffusion terms, see Eq. l |63l l, which are required to resolve 
discontinuities in the magnetic field. However the point is that these 
terms are explicitly added to the SPH calculation and can be turned 
off where they are not necessary (for example using the switches 
described in i)2,2,3| and £|2.3.2t whereas the diffusion present in a 
grid-based code is intrinsic and always present. 



3.3.3 2D: Orszag-Tang test 



The evolution of the compressible Orszag-Tang vortex system 
dOrszag & T ang 1979) involves the interaction of several shock 
waves traveling at different speeds. Originally studied in the 
context of incompressible MHD tur bulence, it has late r been 
extended to the compre ssible case dDahlburg & Picond 1 19891; 
IPicone & Dahlburg 199]]). It is generally considered a good test to 
validate the robustness of n umerical MHD schemes and has been 
used by many authors (e.g.lRvu & Jones|[T995l; iDai & Woodward! 
1 19981 ; Ijiang&Wul 19991: ll.ondrillo & Del Zannall2000l). In the SPH 
context, this test has been discussed in detail by Pried J2004I) and 
IPrice & Monaghanl J2005T) . 

The problem is two dimensional with periodic boundary conditions 
on the domain [0, 1] x [0, 1]. The setup consists of an initially 
uniform state perturbed by periodic vortices in the velocity field, 
which, combined with a doubly periodic field geometry, results in 
a complex interaction between the shocks and the magnetic field. 
The velocity field is given by v — vo[— sin (2ny), sin (27ra;)] 
where vq — 1. The magnetic field is given by B = 
Bo[— sin (27n/), sin (47ra;)] where Bo = 1/V47T. Using the Euler 
potentials this corresponds to a = A z = £>o/(27r)[cos (27ry) + 
| cos (47rx)]. The flow has an initial average Mach number of 
unity, a ratio of magnetic to thermal pressure of 10/3 and we use 
a polytropic exponent 7 = 5/3. The initial gas state is therefore 
P = 5/3Bq = 5/(12tt) and p = jP/v = 25/(36tt). Note that 
the choice of length and time scales differs slightly between vari- 
ous im plementations in the liter ature. The setup used above fo llows 
that of Irvu & Joiiesl dl995T) and lLondrillo & Del Zannal d2000l) . 
We compute the problem using 512 x 590 particles initially placed 
on a uniform, close-packed lattice. The density at t = 0.5 is 
shown in Fig.[l4]using both the standard SPMHD formalism (left), 
see i)2.3.1| and the Euler potential formalism (right), see ^2.3.3 1 
The Euler potential formulation is clearly superior to the standard 
SPMHD method. This is largely a result of the relative require- 
ments for artificial resistivity in each case. In the standard SPMHD 
method the application of artificial resistivity is crucial for this 
problem (that is, in the absence of artificial resistivity the density 
and magnetic field distributions are significantly in error). Using 
the Euler potentials we find that the solution can be computed using 
zero artificial resistivity, relying only on the "implicit smoothing" 
present in the computation of the magnetic field using SPH opera- 
tors in Eqs. ([59} and d60t . This means that topological features in 
the magnetic field are much better preserved, which is reflected in 
the density distribution. For example the filament near the center of 
the figure is well resolved using the Euler potentials but completely 
washed out by the artificial resistivity in the standard SPMHD for- 
malism. Also the high density features near the top and bottom of 
the figure (coincident to a reversal in the magnetic field) are much 
better resolved using the Euler potentials. 

A further advantage of using the Euler potentials is that the field 
lines can be plotted directly as equipotential surfaces of the poten- 
tials. The field lines corresponding to Fig. [14] are thus shown in 
Fig.Ql 

In order to enable a comparison between different codes, we also 
show the ID pressure distribution along the lines y = 0.4277 
and y = 0.3125 in Fig. [131 which may be compared t o simi- 
lar plots given in F ig. 1 1 of lLondrillo & Del Zannal {2000) and in 
|jiang&WuUl999l) . 
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Figure 14. Density distribution in the two dimensional Orzsag-Tang vortex problem at t = 0.5. The initial vortices in the velocity field combined with a 
doubly periodic field geometry lead to a complex interaction between propagating shocks and the magnetic field. Results are shown using 512 X 590 particles 
using a standard SPMHD formalism (left) and using the Euler potentials (right). The reduced artificial resistivity required in the Euler potential formalism 
leads to a much improved effective resolution. The plot may be compared to results shown at comparable resolution in Fig. 14 of Dai & Woodward (1998). 





Figure 15. Pressure distribution in the Orzsag-Tang vortex problem (as 
in Fig. 114) along a ID cut taken at y = 0.4277 (upper panel) and 
y = 0.3125 (lower panel). The plo ts can be compare d to Fig. 11 in 
lLondrillo & Del Zann j fcOQOl) . see also ljiang &Wull 19991) . 



Figure 16. Magnetic field lines (contours of the Euler potential a) in the 
two dimensional Orszag-Tang vortex problem at t = 0.5 (correspond- 
ing to the right panel o f Fig. 114) . Th is may be compared to Fig. 15 in 
iDai & Woodward! £l998) and Fig. 10 in Londril lo & Del Zannd feOOfJ) . 



3. 3. 4 3D: MHD blast wave 

There appear to be very few three dimensional MHD solutions pub- 
lished in the literature. Here we perform an MHD version of the 
Sedov test, identical to the hydrodynamic test with the addition of 
a uniform magnetic field in the x— direction, th at is B = \Bq, 0, 0] 
with Bo = 3.0. A similar test has been used bv lBalsaralfcoOll) for 
testing a 3D Adaptive Mesh Refinement code although with weak 



magnetic fields. Here we perform the test in the strong field regime 
such that the geometry of the blast is significantly constrained by 
the magnetic field, testing both the magnetic field evolution and the 
formulation of magnetic forces in the code. Initially the surround- 
ing material has zero thermal pressure, meaning that the plasma 
/3 P ias is zero (ie. magnetic pressure infinitely strong compared to 
thermal pressure). However, this choice of field strength gives a 
mean plasma /3 p i as in the post-shock material of /3 p i as ~ 1.3, such 
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Figure 17. Results of the MHD blast wave test at t = 0.05 at a resolution of 1 million (100 3 ) particles. Plots show (left to right, top to bottom) density, 
pressure, magnitude of velocity and magnetic field strength (with overlaid field lines), plotted in a cross-section slice through the 2 = plane. 



that the magnetic pressure plays an equal or dominant role in the 
evolution of the shock. 

The initial Euler potentials for the blast wave are: 



-Bqz and = y 



(69) 



where an offset is applied to each potential at the boundaries to en- 
sure periodicity. 

The results of this problem at t = 0.05 are shown in Pig. |17l where 
plots show density, pressure, magnitude of velocity and magnetic 
field strength in a cross section slice taken at z — 0. In addition the 
magnetic field lines are plotted on the magnetic field strength plot. 
In this strong-field regime, the magnetic field lines are not signifi- 
cantly bent by the propagating blast wave but rather strongly con- 
strain the blast wave into an oblate spheroidal shape. The density 
(and likewise pressure) enhancement in the shock is significantly 
reduced in the y— direction (left and top right panels) due to the 
additional pressure provided by the magnetic field which is com- 
pressed in this direction (bottom right panel). 



4 SUMMARY AND OUTLOOK 



drodynamics are derived self-consistently from a Lagrangian and 
account in particular for the so-called "grad-h"-terms. Contrary to 
other approaches, we also account for the extra terms in the gravita- 
tional acceleration terms that stem from changes in the smoothing 
length. This part of the code has been extensively tested on a large 
set of standard test problems. The code performs very well, in par- 
ticular its conservation properties are excellent. While the "grad- 
h"-terms slightly improve the accuracy, in typical applications in- 
volving neutron stars the differences to the older set of equations 
are very minor. 

We evolve the magnetic fields with so-called Euler potentials which 
are advected on the SPH-particles. They correspond to a formula- 
tion of the magnetic field in terms of a vector potential, therefore, 
the V • B = constraint is satisfied by construction. To handle 
strong shocks artificial dissipative terms were introduced in these 
potentials, but for several tests no artificial dissipation is required 
and the corresponding terms can be switched off. The Euler poten- 
tial approach shows in all tests a considerably higher accuracy than 
previous magnetic SPH formulations and is our method of choice 
for our future astrophysical applications of the MAGMA code. 



We have introduced a new, 3D code, MAGMA, for astrophysical 
magnetohydrodynamic problems that is based on the smoothed par- 
ticle hydrodynamics method. The equations of self-gravitating hy- 
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